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Abstract. Reaction-diffusion equations deliver a versatile tool for the description 
of reactions in inhomogeneous systems under the assumption that the characteristic 
reaction scales and the scales of the inhomogeneities in the reactant concentrations 
separate. In the present work, we discuss the possibilities of a generalization of reaction- 
diffusion equations to the case of anomalous diffusion described by continuous-time 
random walks with decoupled step length and waiting time probability densities, the 
first being Gaussian or Levy, the second one being an exponential or a power-law 
lacking the first moment. We consider a special case of an irreversible or reversible 
A —¥ B conversion and show that only in the Markovian case of an exponential waiting 
time distribution the diffusion- and the reaction-term can be decoupled. In all other 
cases, the properties of the reaction affect the transport operator, so that the form of 
the corresponding reaction-anomalous diffusion equations does not closely follow the 
form of the usual reaction-diffusion equations. 
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1. Introduction 

Many phenomena in systems out of equilibrium can be described using a picture of 
react ion- diffusion. Examples can be found in various disciplines, above all in chemistry 
but also in physics, ecology and others. Examples from physics are the trapping and 
annihilation of excitons and the electron-hole recombination in solids. In ecology, 
there are e.g. the predator-pray relations. Both reaction-diffusion with normal and 
anomalous diffusion have been extensively studied over the past decades. However, for 
the latter, a general theoretical framework is still absent. In this article, we discuss 
a special case of the monomolecular conversion under subdiffusion and show that the 
mesoscopic approach leads to equations different in form from what could be regarded 
as a straightforward generalization of the reaction-diffusion scheme. 

The mesoscopic approach leading to reaction-diffusion equations is valid if there is 
a strong scale separation between the typical reaction scale and the size of the system's 
inhomogeneities. The corresponding reaction-diffusion equations (for normal diffusion) 
typically have the form 

= KAQ ± k&?C? ■ ■ ■ cr , (1) 

which simply follows by adding a diffusion term to a classical kinetic equation for 
the corresponding reaction. Here, denotes the diffusivity of the component i, the 
integer powers rij correspond to the stoichiometry of the reaction, and denotes the 
corresponding reaction rate. 

However, many physical systems exhibit anomalous diffusion (see e.g. 121 E] 
for reviews and popular accounts), which is not adequately described by Fick's law. 
Many cases of subdiffusion are successfully modeled within the continuous-time random 
walk framework (CTRWs) with power-law on-site waiting time distributions lacking 
the first moment. These distributions typically have the form w(t) oc with 
< a < 1. Examples include, among others, dispersive charge transport in disordered 
semiconductors, contaminant transport by underground water and motion of proteins 
through cell membranes. On the other hand, successful search strategies in animal 
motion can be described by Levy walks or flights, often in combination with broad 
waiting time distributions. Levy flights are also used as a model for the transport on 
annealed polymer chains [HE], which may be relevant for the gene expression 

For anomalous diffusion, the Fickian diffusion equation is changed for an anomalous 
diffusion equation involving fractional derivatives. For subdiffusion, the equation for the 
concentration C(x,t) of diffusing particles reads 



dC(x,t) 

— M - = K a0 D t 



K a0 Dl- a AC(x,t), (2) 



with the corresponding (anomalous) diffusion coefficient K a , where qD^ stands for the 
operator of a fractional Riemann-Liouville derivative, 
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with n — [/?] + 1 ([x] stands for the whole part of the number x) and v = n — (3. For 
a Levy flight, i.e. the random walk process with the power-law distribution of the step 
lengths, X(x) oc the corresponding equation reads 

^^l = K,A^C(x,t), (4) 

where A M ^ 2 stands for the Riesz symmetric fractional derivative acting on the spatial 
variable 7]. For a "sufficiently well-behaved" function f(x) it can be expressed through 
the Liouville - Weyl derivative [H|: 

for fi ytz 1/2, and for fi = 1/2 through the derivative of the Hilbert transform of /: 

A v* /W = 41 r «m. (6 ) 

dxix i_ 00 x-i 

Reactions under anomalous diffusion were discussed by several authors. However, 
most attention was payed to the description of the elementary act of reaction on the 
microscopic scale (HI UHl HH E21 HHj • Mesoscopic approaches were used e.g. in [llj for 
subdiffusion, where equations of the type 

= K w oA 1-0 * ACi(r, t) + fi (7) 

were postulated for different components in a multi-component system, and in Refs. 
|15| ITo] . where front propagation was discussed for symmetric and asymmetric Levy 
flights, respectively, see also Ref.fTTl and Ref. 0, where a Levy diffusion term was 
added to a "normal" reaction-diffusion equation to describe target search processes on 
the DNA. 

In what follows, we discuss the derivation of the react ion- anomalous diffusion 
equations for a special case of the simple monomolecular conversion A — > B under 
a CTRW transport mechanism (where our approach however differs from the one of 
our previous publication Ref.[18j). We consider subdiffusion, Levy flights and the 
combination of both. Moreover, a reversible conversion A ^ B is also considered. As 
we proceed to show, the Markovian situation of a (symmetric) Levy flight is described 
correctly by the reaction-superdiffusion equation 

^^ = K^ H/2 C t (r,t) + f l (8) 

with Cj being A or B and the reaction terms fi = ±kA. On the other hand, the situation 
for the non-Markovian subdiffusive transport is much more involved. The irreversible 
reaction can be described by an equation for A with the transport term depending on 
the reaction rate, and the equation for the reversible case cannot be casted in a form of 
something resembling a reaction-diffusion equation. 

The article is structured as follows: In Section |21 we derive the equation for the time 
evolution of the educt concentration A in an irreversible reaction. The behavior of the 
product concentration B is discussed in section |3] Section |U is devoted to a mesoscopic 
approach to reversible conversions. 
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2. The educt concentration in the irreversible conversion A — > B 

In what follows, we consider the situation where ^4-particles are converted into B at a 
constant conversion rate k independent on their position. Thus, the survival probability 
of a single A-particle in the time interval [t',t] is — &A(t — t') = exp[— «(£ — £')]. 

We will use one-dimensional notation in the following, the generalization to higher 
dimensions is straightforward. An example for this situation is the decay of a radioactive 
isotope in the groundwater, where the reaction and the transport mechanism are fully 
decoupled. We are interested in the mathematical description of the situation, where 
the transport is given by a decoupled CTRW process with given step length and waiting 
time distribution. Our derivation of reaction-anomalous diffusion equations is parallel 
to the derivation of the pure anomalous diffusion equations in jT]. 

We can put down an equation for the probability density function (pdf) of the 
positions x of the particles, which have just made a jump at time t: 

/OO ft 
I 7] A (x', t')e- K{t - tr) i){x -x',t- t')dx'dt' + A(x, 0)<S(t).(9) 
■oo JO 

Here, 4>(x,t) is the jump pdf given by the probability density in space and time to 
make a jump of length x at time t after the last jump. The meaning of the equation 
is that for whatever t > an A- walker that has just arrived at x has come there 
from some other site, where it had survived as A during the whole waiting time. The 
second term corresponds to the initial condition that at time t = all particles are 
assigned a new waiting time. Here, we have additionally assumed that the jump length 
distribution does not depend on the position of the walker and that the waiting time pdf 
is constant in time and space. Furthermore, ip(x — x',t — t') is assumed to be decoupled 
ij)(x -x',t- if) = X(x - x')w(t - tf). 

In order to get the equation of motion for the A-particles, i.e. for the concentration 
A(x,t), we connect it to T]a(x, t) over 

A(x,t)= [ dt / 7Mfot / )e"* (t "°*(t-0: (10) 
Jo 

where ty(t — t') is the probability to stay at site x for a time (t — t') after a jump. It is 
given by 

V(t) = 1 - I dt'wit'). (11) 
Jo 

Both Eqs. (JHJ) and contain convolution integrals and can be solved by Fourier-Laplace 
transform. Using the shift theorem for the Laplace transform, we get 

!(*,«)= . (12) 

(u + k) [1 — lf)(k, u + k)] 
Before we can return to the space- and time-domain, we have to specify the jump 
pdfs. We are interested in the continuum limit of the equations corresponding to large 
scales and long times, i.e. to (k,u) — > (0,0). A characteristic function of a Gaussian 
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jump length pdf with variance 2a 2 will then be approximated by X(k) ~ 1 — k 2 a 2 . 
A characteristic function of a broad Levy distribution, X(k) = exp(— <7 M |fc| M ), can be 
approximated through A(A;) — 1 — a^\k\^. For a broad waiting time pdf of a Pareto 
(power-law) type, w(t) ~ ar a t~ 1 ~ a , one infers the following asymptotics in Laplace 
space using a Tauberian theorem, w(u) ~ 1 — T(l — a)u a T a . For the Markovian case, 
as exemplified by the exponential waiting time pdf, w(t) = r -1 exp(— t/r), one has 
w(u) ~ 1 — ut in the continuum limit, which corresponds to a = 1. Eq. (jl2J) then can 
be rewritten in the following form 

uA(k, u) - A(k, 0) = -{u + K f- a — lifer A(k, u) - nA(k, u) (13) 

simplifying the inverse transforms. For the inverse Fourier transformation, we use 
J r ~ 1 {-k 2 f(k)} = Af(x), and ^{-{k^ f(k)} = A^ 2 f(x). Moreover, we introduce 
the notation K^ a = cr^[r a r(l — a)] -1 for what later will be identified as the generalized 
diffusion coefficient. The inverse Laplace transform of the left hand side of the equation 
is simply the first time derivative, since C~ l {ug{u) — g(0)} = dg{t)/dt. 

We first combine the Gaussian and Levy distributed jump length pdf with an 
exponential waiting time pdf. In this case, the pre-factor of A(k, u) in the first term 
on the right side of the equation does not depend on u. After inverse transforming 
the equation, it becomes a time-independent operator acting on the concentration as a 
function of the coordinates. For a Gaussian jump length distribution, our equation (fH?j) 
now reads 

dA(x,t) 



dt 

and for Levy flights, 
dA(x,t) 



K 2 ,iAA(x, t) - kA(x } t), (14) 



K^A^ 2 A(x, t) - kA(x, t). (15) 



dt 

Hence, the separation of the transport- and the reaction-term is perfectly exact. 
For Levy flights, the Laplace operator is just changed for the Riesz-Weyl fractional 
derivative. 

Now, we turn to subdiffusion and consider a Gaussian distribution of the step 
lengths (fi = 2) combined with a broad waiting time pdf of a Pareto type with < a < 1. 
From Eq.([13|) we then get 

dA ft^ = K 2 , a T t {l - a, k)AA(x, t) - kA(x, t), (16) 
with the transport operator %(1 — a, k)A, which is now time-dependent, 



1 fd e~ K (*-*') 



Its form follows from the shift theorem for the Laplace transform. We see that the 
reaction parameter enters the transport-term, and the transport operator %(1 — a,K) 
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reduces to a fractional derivative only for k = 0. Using the Laplace transform property 
of the Riemann-Liouville fractional derivative, C~ 1 {u~ a f(u)} = oD 1 T a f(t) for a > 0, 
and using the shift theorem, the temporal part of a transport operator (in its action on 
the arbitrary function of time f(t)) can be transformed into a form JH] 

T t (l - a, K)f(t) = exp(-Kt) D l t - a {exp(Kt)f(t)}. (18) 

One can also easily formulate the equations for the combination of Pareto waiting times 
and Levy jumps being of the form 

= K^ a T t {l - a, k)A^ 2 A(x, t) - kA(x, t), (19) 

with A^/ 2 denoting the symmetrized (Riesz-Weyl) spatial fractional derivative. 

By the way, as shown in j20j > an external force field can be included in the model over 
an asymmetric jump length distribution leading to a fractional Fokker-Planck equation 
with the time fractional operator changed for our operator %(1 — a,K) and with an 
additional reaction term. 



3. Equations for the product concentration 

Let us turn to the equation for the concentration of the S-particles. One can distin- 
guish two different cases: (i) Either a 5-particle takes over the waiting time of the 
A-particle that it was converted from, or (ii) we assign it a new waiting time when it is 
produced. The former means that the conversion is just a relabeling from the standpoint 
of diffusion and that conversion and transport are totally independent. The latter is 
appropriate when A- and U-particles have different diffusive properties, e.g. when they 
are trapped by different kinds of molecules. Then, transport and conversion are partly 
coupled. 



(i) The first case corresponds to the following approach 



r) B (x,t) 



dx' / dt' 



r, B (x',t')+ VA (x',t') (l-e-"(*-*'>) 
x if;(x -x',t- f) \ + B(x, 0)5(t), 



x 



(20) 



which expresses the fact that a 5-particle that has just landed at x at time t has come 
from a site x' at prior time t', where it had either come already as a -B-particle or where 
it had been converted from an A-particle. For the concentration of S-particles, we have 



B(x,t) 



dt' 



V B(x,t') + VA (x,t') (l - e"^-*'))] *(t - 



(21) 



with ^(t) from Eq. ffTTj) . Now, a -B-particle that is at site x at time t has come there at 
a prior time t' either already as a S-particle or as an A-particle and has been converted 
in (t — f). Eq. (}2*Tj) can also be solved using Fourier-Laplace transform and Eqs.©, (fTB|) 




Figure 1. Shown are the the concentrations of A-particles (solid lines) and B- 
particles (dashed lines) for subdiffusion with conversion. The correct results (solution 
of Eas. (|16fl and l|23|l) are shown without dots. They are compared to the solutions 
of the decoupled equations J7J shown with dots. The parameters are: a = 0.75, 
k = 0.001, Kl ~ 7.76 • 10~ 3 . The times shown are t = 200 (left) and t = 2000 (right). 



and pOjl . First, we get 

g( M)+ l( M) = g(M±i(Ml^M, (22) 

which is essentially the Fourier-Laplace transformed subdiffusion equation for the sum 
of the the concentrations C(x,t) = A(x,t) + B(x,t). This is due to the fact that we 
have assumed a complete independence of the transport and the conversion and can 
already be seen adding the two approaches Eqs.tJH} and (J20|) . Using the corresponding 
solutions for the concentration of A-particles, for a Poissonian waiting time pdf, one 
infers an equation of the form ((IJ) or (JSJ). For a power-law waiting time pdf and the 
initial conditions A(x, 0) = 5(x), B(x) = we get 

= K 2>a Dl~ a AB(x, t) + kA(x, t) + 

ot 

+ K 2 , a [ D]- a - %{l - a, «)] AA(x, t). (23) 

The change of the concentration of the i?-particles depends on the concentration of the 
v4-particles at all previous times. This is due to the fact that the I?-particles are already 
"aged" when produced and have a memory for the last jump they have made as an 
A-particle because of the non-Markovian nature of the waiting time pdf. As mentioned 
above, the combination with a Levy distributed jump length pdf leads to the same result 
as Eq. (f2*3*|) with the Laplace operator just changed for its fractional generalization. 

In FigG] we compare the correct solutions, i.e. the solutions of Eqs. (ll6j) and 
(J2HJ), with the solutions of the special cases of Eq.((ZJ) for the conversion. We note 
an even qualitative difference, so the latter justified only by analogy to the normal 
diffusion case cannot be used as an approximation of the exact equations. In order 
to get these results, we did not actually have to solve Eqs. (fTH|l and because we 
could specify the solution from the fact that C(x,t), the sum of the concentrations of 
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A- and -B-particles, fulfills a pure subdiffusion equation. For the conversion reaction 
with the reaction independent on the transport, the concentrations are just given by 
A(x, t) = C(x, t) exp(— K,t), B(x,t) = C(x, t)[l— exp(— k£)], namely by the product of the 
overall particle concentration and the survival probability or the conversion probability, 
respectively. The solution of the pure subdiffusion equation for C(x,t) and the initial 
condition, C(x,0) = 5(x) is known. It is the Fox's H-function, 



C(x,t) 



1.0 



\X\ 



(1 -a/2, a/2) 
(0,1) 



(24) 



The Fox's H-function can be calculated using a series expansion [1. The equations of 
the form Eq.([7]l were solved using a modification of a numerical scheme presented re- 
cently by Yuste et al |2T]. The scheme is a combination of a forward-time-centered-space 
discretization and the Griinwald-Letnikov form of the fractional derivative. 

(ii) Let us now consider the second case and assume that -B-particles are assigned 
a new waiting time at production. Here, we expect to get a decoupled equation of the 
form (j2J) because the past as an A-particle is "forgotten" . We have to start from 



r) B (x,t) 

and 

B(x,t) -- 
This leads first to 

B(k,u) 



oo JO 



r) B {x', t')il){x -x',t- t')dx'dt' + 



+ KA(x,t) + B(x,0)6(t), 
rj B {x,t f )^(t - t')dt'. 
l-w{u) nA{k,u) + B{k,0) 



(25) 
(26) 

(27) 



u l-ip(k,u) 

and then with a Gaussian jump length pdf and the same initial conditions as above to 



dB(x,t) 
dt 



K 2 , aB Dl- aB AB(x,t) + KA(x,t), 



(2f 



the expected decoupled equation. We have denoted the diffusion exponent as in 
order to emphasize that it is possibly different from the one for the A-particles. By 
the way, instead of the reaction-term kA of the conversion we could have an arbitrary 
reaction term that does not depend on the product concentration. 

4. Reversible A ^ B reaction 

Now, we turn to the case of a reversible conversion. We assume that no new waiting 
time is assigned when a particle is converted. We denote the forward reaction rate by 
K\ and the backward rate by Then we have to start from 

,-K 2 (t-t')' 



r} A {x,t) 



oo JO 



X 



x ip(x -x',t- t')dx'dt' \ + A(x, 0)5(t). 



(29) 
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An A-walker that arrives at x at time t has come from another site x' at a prior time 
t', where it had come already as an A-particle and was not converted, or where it had 
come as a B-particle and was converted. For the concentration, we have 

A(x, t)= [ dt'r] A (x, t')e~ K ^ t,] m{t - t') + 



+ dt'r] B (x, t') (l - e -^(*-*')) q?(t - t'). (30) 

An A-particle at site x at time t has come to this site already as an A at time t' and 
has not been converted and moved since, or it has come there as a .B-particle, was 
converted and has not moved in the mean-time. Because of the "symmetry" of the 
reaction, the equations for the 5-particles can be directly inferred from the equations 
for the A-particles. We can still perform Fourier-Laplace transform. Using a Gaussian 
jump length pdf, an inverse power-law waiting time pdf, and the initial conditions 
A(x, 0) = <5(x), B(x, 0) = 0, we find for the A-particles 

K 2 , a k 2 + [(u + «i) - (u + kx) 1 -^ ] x 



A(u, k) 



[K 2 , a k 2 + (u + K2) a ][K 2 , a k 2 (u + Ki y- a + (u + Kl )] + 



X [U a - 1 - (U + K2)*- 1 ] + (U + K 2 ) 



(31) 

+ [u a (u + Ki) 1 " - [U + K X )][{U + K 2 ) a - U a ] ' 

However, after the inverse Fourier-Laplace transform, the equation of motion does not 
take any simple form, let alone the form of a react ion- diffusion equation. The decoupled 
scheme, Eq.©, corresponds to a different equation, 

A( M = K ^ k2 + K ^ + u ° (v) 

{U ' ' [K 2 ^ + u a }[K 2 , a k 2 u^ + u + Kl + K2 y 1 ' 

We have made some simple simulations in order to test how the decoupled equations 
perform for this case. For the conversion, the random walkers are independent, and we 
can simply repeat the random walk procedure many times (10 6 times). We used the 
power law waiting time pdf with a cut-off at small times guaranteeing the normalization, 
w(t) = ar a t~ 1 ~ a for t > r and w(t) = otherwise. The conversion is independent of 
the jumps and takes place with a constant probability Pa,(b) = [1 — exp(— Ki^At)} in 
each time-step of length At. We can see in FigO that the coincidence with the correct 
result is somewhat better for large times than in the case of an irreversible conversion. 



5. Conclusions 



We discussed generalizations of the reaction-diffusion scheme to the case of anomalous 
diffusion for a special case of a simple conversion reaction A — > B or A ^ B. Although 
the reaction and the particle transport were assumed to be independent, the reaction- 
term and the transport-term do not separate in the case of subdiffusion. This means 
that the transport operator in the corresponding equations depends on the properties 
of the reaction. The simple equations with separated reaction- and diffusion-terms are 
not exact. Comparing the exact solution with the equations with decoupled reaction- 
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Figure 2. Shown is the result of the simulation: A-particles (solid line) and B- 
particles (dashed line), and the numerical solution of the decoupled equations, Eq.Q 
or Ea. 1(3 2(1 : A-particles (squares) and B-particles (dots). The parameters are: a — 0.75, 
K a = 0.0138, ki = 0.01, n 2 = 0.001. The times shown are t = 200 (left) and t = 2000 
(right). 

and diffusion-terms shows that the latter deliver a rather poor approximation for the 
case of an irreversible reaction and perform somewhat better in the reversible case. 
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